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ABSTRACT 

The gas equation of state may be one of the critical factors for the disk insta- 
bility theory of gas giant planet formation. This letter addresses the treatment of 
H 2 in hydrodynamical simulations of gravitationally unstable disks. In our dis- 
cussion, we point out possible consequences of erroneous specific internal energy 
relations, approximate specific internal energy relations with discontinuities, and 
assumptions of constant L\. In addition, we consider whether the ortho/para ra- 
tio for H 2 in protoplanetary disks should be treated dynamically as if the species 
are in equilibrium. Preliminary simulations indicate that the correct treatment 
is particularly critical for the study of gravitational instability when T = 30-50 
K. 



Subject headings: equation of state - hydrodynamics 
processes - planetary systems: protoplanetary disks 



1. INTRODUCTION 



instabilities - molecular 



Researchers disagree on several key issues of disk evolution, particularly regarding frag- 
mentation and planet formation, despite numerous studies of gravitational instabilities in 
protoplanetary disks (see Durisen et al. 2006 for a review). Recent three-dimensional radia- 
tive hydrodynamics simulations of protoplanetary disks with gravitational instabilities (GIs) 
demonstrate disparate evolutions, and the differences involve the importance of convection 



20061 ). 



(Boss 


2005; 


Cai et al. 


2006; 


Bolev et al. 


2006; 


Maver et al. 
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These disparities may be due to differences in the treatments of radiative transfer and 
the assumed equation of state. Our own group has undertaken systematic improvements 
to the Indiana University Hydrodynamics Group code in these two areas. Boley et al. (in 
preparation) will present a new scheme for 3D radiative cooling in disks that passes a series 
of accuracy and reliability tests, which are suggested as a new standard for disk stability 
studies. Because the equa tion of state can severely change the outcome of a simulation 



( iPickett et al.lll998l . l2000al ). we have also implemented an internal energy that takes into 
account the translational, rotational, and vibrational states of H 2 . During its development, 
we noticed that in all treatments to date of planet formation by disk instability, the effects 
of the rotational states of H 2 have been, at best, only poorly approximated. The purpose of 
this Letter is to draw attention to possible consequences of various approximations for the 
internal energy of H 2 that are in the literature and to disseminate the correct treatment for 
H 2 . We have begun a series of simulations to investigate the severity of using the incorrect 
internal energies. We mention only preliminary results in this Letter, but the simulations 
will be fully discussed in a future paper. 

This Letter is laid out as follows: In §2 we outline the correct treatment of H 2 under 
simple conditions, i.e., no dissociation and no ionization. We discuss the ortho/para hydrogen 
ratio for optically thick protoplanetary disks in §3, and we summarize the main points of 
the letter in §4. 



2. THERMODYNAMICS 



In this section we su mmarize relevan t thermodynamic relations. For an in-depth discus- 
sion, we refer the reader to lPathrial (119961 ). Consider the following thermodynamic properties 
of an ideal gas: Let E be the internal energy for N particles, e the specific internal energy, e 
the internal energy density, p the pressure, T the gas temperature, p the gas density, \i the 
mean molecular weight in proton masses, c v the specific heat capacity at constant volume, 
Z the partition function for the ensemble, z the partition function for a single particle, and 
R = k/m p , where k is Boltzmann's constant and m p is the proton mass. We only consider 
independent contributions to the partition function from translation, rotation, and vibration 
represented by Z = Z tTain Z mt Z v fo = z N = (^tran^rot-Svib)^- The internal energy E and the 
specific internal energy e can be calculated by 

,<9 In z 



E = NkT 



dT 



R T 2 dlnz 
H dT 



(1) 



for constant p. Because the gas is ideal, c v = de/dT. If c v is constant, then e = c v T. 
Black fc Bodenheimerl (119751 ) calculate c v from the Helmholtz free energy, which is valid, 
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but then assume e = c„T, which is invalid, bec ause c v is dependent on T. Other authors 



have followed suit (e.g., Whitehouse fc Bate 20061 ). and this assumption could be troublesome 



for gas dynamics in a hydrodynamics simulation (see below). We do note that the severity 
of this error may depend on the state variables evolved in a given code, and only the authors 
who employ e = c v T will be able to say in detail how it affects their simulations. 



2.1. Molecular Hydrogen 



Molecular hydrogen exists as para hydrogen and as orthohydrogen where the proton 
spins are antiparallel and parallel, respectively. The partition function for parahydrogen is 



^(2j + l)exp(-j(j + l)0 rot /T) 



(2) 



and the partition function for orthohydrogen is 



^3(2j + l)exp(-j(j + l)# rot /T) 



(3) 



Jodd 



where 8 mt = 85.4 K (IBlack fc Bodenheimerlll975l ). When the two species are in equilibrium, 
z p + z a . However, the ortho/para ratio (b:a) could also be frozen if no efficient mech- 



^rot 



(a/(a+b)) r(b/(a+b)) 



anism for converting between the species is available. This leads to z TOt = z v ' ' " z 
where z' Q = z Q exp (26 mt /T). The additional exponential is required in the orthohydrogen 
partition function when the ortho and para species are at some fixed ratio to ensure that 
rotation only contributes to the internal energy once the rotational states are excited, i.e., 
z' — ► 1 as T — > 0. 

To consider the vibrational states, we approximate the molecule as an infinitely deep 
harmonic oscillator, where 

1 - exp (-0 v ib/T) 



#vib = 5987 K (jDraine et al.lll983l ). Because we are only interested in temperatures T < 1500 
K where dissociation of H2 is insignificant, we can ignore differences between equation (4) 
and a proper z v ^, which would take into account the anharmonicity of the molecule and that 
the molecule has a finite number of vibrationally excited states. 



We can use equation (1) to write the specific internal energy for H 2 



(H 2 ) 



R 
~2 



3 r r T 2 dz mt exp 
2 T + l~ t ^T + 9vih 



7 vib 



1 - exp (-0vib/T) 



(5) 
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When the gas is ideal and dissocia tion and ionization can be ignored, the first adiabatic 
exponent Ti — 1 + R//j,c v = c p /c v = 7 (jCox fc Giuh1ll968l ). Figure 1 indicates the r x profiles 
fo r the equilibrium case, pure parahydrogen, and a 3:1 mix, which is consistent with Figure 2 
of iDecampli et al.l (119781 ). as it should be because our derivation of c v is equivalent to theirs. 
The specific internal energy profiles as calculated by equation (5) are shown in Figure 2 by 
solid lines. Curves for c v T are also shown in Figure 2 by dashed lines. The approximation 
e = c v T gives quite different behavior from the correct e, e.g., the incorrect 3:1 curve most 
closely follows the correct pure parahydrogen curve. The offset of the correct 3:1 mix profile 
is due to the energy stored in the parallel spins of the protons. The dotted line represents 
the curve e = 3/ 4RT f o r T < 100 K and e = 5 /4RT for T > 100 K, which is the energy 



equation used by iBossI (119841 . l200ll . 12002k 120051 ). Although this may be a reasonable ap- 
proximation for e when the hydrogen species are in equilibrium, this approximation could 
be troublesome for the hydrodynamics. This is because the discontinuity guarantees that 
c v becomes very large near 100 K, which forces T\ to unity. Such discontinuities may ar- 
tificially favor gas gia nt formatio n by d i sk instability, becau se isothermal (Ti) gas behavior 
favors frag mentat i on (|Bossl 119971 . |2000| ; iPickett et al.lll998l ). Indeed, the clumping that is 

y, a constant Ti 



2005 



reported by Boss! (2001) occurs at T m 100 K (see Boss's Figure 2). Final 



approximation ( Pickett et al. 20031: Lodato &: Ricel 2004 - L Mayer et al. 



Mayer et al. 



200 



6: 



Mejfa et al. 



200 



5: 



Cai et al. 



200C 



2004 



Rice et al 



ft iBoley et aPhood : iMaver et al. 



2003, 



20061 ) poorly represents e in the temperature regime where Jupiter probably formed; neither 
ri = 5/3 or 7/5 can be assumed confidently. 

The dynamical effects that could result from assuming e = c v T are evaluated by taking 
the temperature derivative of the dashed curves in Figure 2. The results are displayed in 
Figure 3, and the curves depart drastically from the correct Ti profiles . As mentioned above, 



c„ can be calculated from the Helmholtz free energy as is done by iBlack fc Bodenheimer 



(119751 ) . This makes it possible to compute Ft from c v correctly but then evolve the gas 



(e.g., 


B 


ack & Bodenheimer 


1975: 


Pickett 


1995|; 


Wadslev et al. 


2004) 



Bosslll984j . l200ll : lMonaghanlll992l : IStone fc Normanlll992 



the ideal gas law p — Rf pipT is assumed as well, then 
effective Ti profiles like those shown in Figure 3 seem to be unavoidable when e = c v T is 
assumed for a temperature dependent c v . Because fragmentation becomes more likely as 
becomes smaller (Rice et al. 2005; Michael et al. in preparation), the e = c v T assumption 
should artificially make fragmentation more likely in some temperature regimes and less 
likely in others. 

Preliminary simulations of a disk with solar composition indicate that when GIs activate 
between 30 and 50K for an equilibrium ortho-para mixture, the e = c v T simulation evolves 
more rapidly and has a more flocculent spiral structure than the correct e simulation for 
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the same cooling rates. In addition, denser substructures form in some spiral arms of the 
e = c v T simulation throughout the simulation, while dense substructures only form during 
the burst of the correct e simulation^. When the instabilities occur outside this temperature 
regime, the differences are diminished. 



3. THE MOLECULAR HYDROGEN ORTHO/PARA RATIO 



As indicated in §2, the dynamical behavior of the gas is dependent on the ortho/para 
ratio and whether the species are in equilibrium for all T. This ratio for various astrophysica 



conditions has b een addressed by several aut h ors fe.g..lOsterbrocklll962l;lDalgarno et al. 



Rodriguez-Fernandez et al 



2000 



Flower et al.l 



Decampli et~ai1 Il978l: [Flower fc Wattl Il984l: [Sternberg fc Neufeldl Il999l: iFuente et al" 



1973 



1999 



20061 ). typically in the context of interstellar 



clouds or photodissociation regions. However, for plausible S olar Nebula condition s, the 
ortho/para ratio has been inadequately addressed; for example, iDecampli et al.l (119781 ) used 
an estimate for the H + number density that was derived originally to give the total gas 
phase ion number density in gas for which dissociative recombination dominates the removal 
of ions. At protoplanetary disk number densities, however, ion removal should be primarily 
on grain surfaces if the ratio of grain surface area to hydrogen nucleon number density is the 
same as it is in diffuse interstellar clouds. 

For protoplanetary disk conditions, the conversion between ortho and parahydrogen 
is principally due to protonated ions such as HJ . Consequently, we will assume that all 
ionizations lead to Hjj~ formation. Another possible conversion mechanism is through inter- 
actions between H 2 and grains; howev er, this conversion might only be significant when the 
temperature drops below about 30 K (ILe Bourlotll2000i ). 



Consider the balance between Hg production by cosmic rays (CR) and Hjj" depletion by 
dust grains: 

(n (H 2 ) = n g Tia 2 vn u (6) 



where n (H 2 ), n i; and n g are the H 2 , ELt , and grain number densities, respectively, ( is the 
ionization rate by CRs and other energetic particles (EP), a is the average radius of the 
grains, and v is the thermal velocity of H^. If we assume standard interstellar extinction, 
then a = n (H 2 ) /n g iia 2 « 10 21 cm -2 , but as we discuss below, this number is ambiguous. 
The ( appropriate for a protoplanetary disk is also ambiguous. Cosmic rays and stellar 



1 The evolution of these simulations can be viewed at http://westworld.astro.indiana.edu/. Click on the 
link titled "H2 ortho-para equilibrium tests" under the "Movies" tab. 
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EPs are important in ionizing the disk surface (jDeschl 12004 ; iDullemond et al.l 120071 ). but 
because these particles are att enuated exponentially with a scale length of about 100 g cm -2 
( jUmebayashi fc Nakanol Il98ll ). stellar EPs probably do not contribute to n^. Moreover, 
protostellar winds could lead to a significant reduction of ( in analogy to CR modulation 
by the solar wind (jWebberlll998l ). However, at a surface density of roug hly 380 g em" 2 , EP 



production by 26 Al decay is as important as CRs with ( ~ 10 19 s 1 (jStepinskil 119921 ). It 



is likely that 10~ 19 s -1 < C < 10~ 17 g~ x . For our estimate we adopt the interstellar rate 
( = 10~ 17 s -1 ( jSpitzer fc Tomaskolll968l ). Using these numbers in equation (6) and adopting 
a thermal velocity of 1 km s -1 , rii w 0.1 cm -3 . By adopting a collisiona l rate coefficient 
a = 1 x 10~ 9 cm 3 s _1 for the interaction with H 2 (jWalmsley et al.l 12004 ). the lower limit 
timescale for ortho and para hydrogen to reach equilibrium is t e = (arii)~ l = 300 yr. 

The equilibrium timescale is short enough that the ortho/para ratio can thermalize in 
the lifetime of a disk, but the equilibrium timescale is longer than the dynamical timescale 
inside about 40 AU: ortho and parahydrogen should be treated as independent species for 
hydrodynamical simulations of young protoplanetary disks. 

What ortho/para ratio should a dynamicist assume for gravitationally unstable proto- 
planetary disk s imulations? The answ er is uncertain. Vertical and radial stirring induced 
by shock bores (IBoley &: Durisenl 120061 ) . which could possibly lead to mixing of the low al- 
titude disk interior with the high-altitude photodissociation region in the disk atmosphere 
(IDullemond et al.l 120071 ). will transport gas through different temperature regimes on dy- 
namic timescales. This could lead to nonthermalized ortho/para ratios like t hose that are 



measured from H_2 rotational transition lines in some photodissociation regions dFuente et al. 



19991 ; iRodrfguez- Fernandez et al.ll2000l ) and in Neptune's stratosphere (IFouchet et al.ll2003l ). 
Moreover, accretion of the outer disk will bring material with a cold history into warmer 
regions of the disk. It is unclear whether the ortho/para ratio of, say, 15 K gas will be ther- 
malized with z /z p » or whether the ortho/para ratio w ill be 3:1, which is the expected 
ratio for H2 formation on cold grains (IFlower et al.ll2006l ). Unfortunately, the ortho/para 
ratio may be critical to the evolution of a protoplanetary disk. As can be seen in Figure 1, 
the pure parahydrogen mix has a Ti that approaches 4/3 for T w 160 K. This could make 
the 160 K regime the most likely region of the disk to fragment because, as Ti decreases, it 
becomes harder for the gas to support itself against local gravitational and hydrodynamic 
stresses (Rice et al. 2005; Michael et al. in preparation). Hydrodynamicists need to consider 
ortho/para ratios between pure parahydrogen and 3:1 because of our ignorance of this ratio 
in protoplanetary disks. 

The above discussion is based on the assumption that o m 10 21 cm -2 , which is probably 
reasonable for very young protoplanetary disks but may not be reasonable for disks with 
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ages of about 1 Ma. Grain growth and dust settling may significantly lower the value of a by 



depleting the total grain area (e.g.. lSano et al.ll2000l ). Because models of T Tauri disks must 



take into acc ount the effects of grain growth in order to m atch observed spectral energy 
distributions (ID'Alessio et ajj [200l| , 120061 ; iFurlan et al.l 120061 ) . there may be a period in a 
disk's evolution when the ortho and parahydrogen change from dynamically independent 
species to species in statistical equilibrium. Suc h a transition may also take place at certain 
radii in a disk, e.g., near edges of a dead zone (lGammidll996l ). As indicated by Figure 1, a 
transition to statistical equilibrium could have significant dynamical consequences for disk 
evolution and may induce clump formation by GIs. 



4. SUMMARY 

The effects of the rotational states of H2 must be explicitly modeled in hydrodynamical 
simulations of protoplanetary disks. Constant Y\ approximations are insufficient for modeling 
the dynamic behavior of the gas because they ignore the transition from a ^ = 5/3 gas to a 
Y\ = 7/5 gas. This transition probably took place near Jupiter's location in the young Solar 
Nebula. Discontinuous e(T) assumptions or the e = c v T approximation could lead to severe 
errors in the behavior of the gas, including artificially driving ri — ► 1. 

Our estimates indicate that ortho and parahydrogen should be treated as separate 
species in hydrodynamics simulations of young protoplanetary disks because the timescale to 
thermalize the ortho/para ratio is longer than a dynamic timescale. However, a reasonable 
approximation to the ortho/para ratio in young protoplanetary disks is ambiguous. More- 
over, as a disk evolves, there may be a transition in the behavior of ortho and parahydrogen 
from independent species to statistical equilibrium, which might trigger disk fragmentation. 
Hydrodynamicists should consider a range of ortho/para ratios and the equilibrium case 
until this ambiguity is resolved. 
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Fig. 1. — Profiles of r\ for an equilibrium mix (black), parahydrogen (red), and a 3:1 or- 
tho/para rat io mix (blue). The dotted line indicates I\ = 4/3; this figure is similar to 
Figure 2 of (IDecampli et al.lll978l ). [Profiles of Ti for an equilibrium mix (solid), parahy- 
drogen (short dash), and a 3:1 ortho/para rat io mix (long dash). T he dotted line indicates 
Ti = 4/3; this figure is similar to Figure 2 of ( IDecampli et al.lll978l ).] 
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Fig. 2. — Specific internal energy profiles. The solid lines indicate the e profiles as calculated 
by equation (5), while the dashed lines indicate c v T. The dotted line indicates a disconti nuous 
e profi le where e = 3/4RT for T < 100 K and e = 5/4RT for T > 100 K, as used by [Boss 
(I1984I ). The inset shows the behavior of the profiles at higher temperatures with the same 
units for the ordinate and abscissa. 
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Fig. 3. — Profiles of F± calculated by taking the temperature derivative of the dashed curves 
in Figure 2. These profiles vary largely from the correct profiles in Figure 1. The equilibrium 
curve's maximum is at about Fx = 2.3. 



